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ABSTRACT 

We derive an upper limit of 3 mJy (95% confidence) for the flux density 
at 317 MHz of the Geminga pulsar (J0633+1746). Our results are based on 
7 hours of fast-sampled VLA data, which we averaged synchronously with 
the pulse period using a period model based on CGRO/EGRET gamma-ray 
data. Our limit accounts for the fact that this pulsar is most likely subject 
to interstellar scintillations on a timescale much shorter than our observing 
span. Our Bayesian method is quite general and can be applied to calculate the 
fluxes of other scintillated sources. We also present a Bayesian technique for 
calculating the flux in a pulsed signal of unknown width and phase. Comparing 
our upper limit of 3 mJy with the quoted flux density of Geminga at 102 MHz, 
we calculate a lower limit to its spectral index of a ~ 2.7 (F(u) oc v~ a ). We 
discuss some possible reasons for Geminga's weakness at radio wavelengths, and 
the likelihood that many of the unidentified EGRET sources are also radio-quiet 
or radio-weak Geminga-like pulsars. 



1. Introduction 

Geminga has long been known as a strong X-ray and gamma-ray continuum source 
( [Fichtel et al. 1975|) and, in 1987, was identified with an optical star (pignami et al. 1987| ). 
In 1992, Geminga was detected at X-ray ( [Halpern fc Holt 1992] ) and gamma-ray ( |Bertsch| 



et al. 1992|) energies as a 237-ms pulsar with a period derivative of 10.95 x 10 15 s s" 



indicating a characteristic age of 3.4 x 10 5 yr and a magnetic field of 1.6 x 10 12 Gauss. A 
proper motion of 170 ± 6 mas yr^ 1 and a parallax of 6.36 ± 1.7 mas have been measured for 
Geminga's optical counterpart ( |Caraveo et al. 1996|) . From these measurements, a distance 



of 157±|1 pc and a transverse velocity of 140 ± 15 km s 1 have been calculated. 
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Since its discovery as a gamma-ray source in 1975, many observers have unsuccessfully 
attempted radio detection of Geminga both as a continuum source flSpelstra fc Hermsen| 
|1984[ ) and as a pulsar ( jSeiradakis 1992| ). Recently, however, three independent groups 



( iMaloleev & Malov 1997| ; |Kuz'min & hosovskn lW7| ; fthitov & Pugachev 1991 ) have 



reported successful detections of pulsed radio emission from Geminga. All 3 detections 
were made at the Pushchino Radio Astronomy Observatory at a radio frequency of 102 
MHz, a lower frequency than those used for previous searches. All three groups calculate a 
dispersion measure (DM) of roughly 3 pc cm -3 . The measured fluxes vary greatly, perhaps 
due to interstellar scintillations, and range from 5-500 mJy Reported pulse widths range 
from 10 to 120 ms. Recently, Shitov & Malofeev (1997) reported a 41 MHz detection 
of Geminga with mean flux density 300 mJy. Comparing this result with Malofeev and 
Malov's quoted 102 MHz mean flux of 60 mJy, we may estimate a spectral index of a ~ 1.8, 
where F{y) oc u~ a . Comparing Malofeev and Malov's quoted mean flux of 60 mJy with the 
previous 1 mJy upper limit at 1.4 GHz ( |Seiradakis 1992| ), we calculate a lower bound on 
Geminga's spectral index of a ~ 1.6. 

Since Geminga is the only known gamma-ray pulsar which is not also a strong radio 
source, confirming its existence as a radio pulsar and determining the reasons for its 
radio weakness are essential to understanding the relationship between pulsar radio and 
gamma-ray emission. Geminga is also important as it may be a prototype for a class 
of radio-weak or radio-quiet high-energy pulsars which could account for some of the 80 
EGRET sources not yet identified with known pulsars or active galactic nuclei. 

We have therefore revisited the search for Geminga as a radio pulsar with a long, low 
frequency, fast-sampled VLA observation. In this paper, we present our results, introducing 
Bayesian methods for determining the pulsed flux in a signal of unknown width and 
phase and for estimating the intrinsic flux density of a source, accounting for interstellar 
scintillations. 



2. Data 

We observed Geminga (PSR J0633+1746) on 2 February 1998 with the D-configured 
VLA in phased array mode. The data spanned a bandwidth of 25 MHz, centered on 
317.5 MHz, with 14 (not necessarily consecutive) 1-MHz channels (situated to avoid 
known interference frequencies). Full polarization data were recorded with the High Time 
Resolution Processor ( Moffett 19971) using a sampling frequency of 1152 Hz. We obtained 
7 one-hour observations (scans) of Geminga which spanned a total time of 9 hours. Two 
other known pulsars (B0329+54 and B0950+08) were observed for 10 minutes each with 
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the same setup as that used for Geminga. Each scan was started on a 10 second tick to 
ensure accurate pulse phase referencing between scans. This 10-second tick was tied to the 
VLA's hydrogen maser, which is compared to Universal Time through the GPS. The array 
was phased, with rms phase errors less than 1°, between consecutive scans. We observed 
a strong VLA calibrator source (B0813+482) with known flux density (S = 45 Jy) to 
calculate the gain in each of the 56 channels (14 frequencies and 4 polarizations). As we 
expect only a 7-mJy uncertainty in our flux density calibration due to radiometer noise, 
our main source of uncertainty is the intrinsic source flux variation. By comparing the 
VLA flux measurement with those from the Westerbork and Texas surveys ([Rengelink et 



al. 1997| ; |Douglas et al. 1996Q , we estimate a maximum intrinsic flux density variation of 
10%. We therefore expect errors in our final reported flux densities to be less than 10%. 



3. Analysis 

For each of the 7 Geminga data sets, we summed polarizations to calculate the 
total intensity and dedispersed all 14 frequency channels using a DM of 3.2 pc cm -3 , the 
midpoint of the DM range predicted by the TC ( |Taylor fc Gordes 1993| ) model for Galactic 



electron density given Geminga's measured distance and direction. Our estimated DM is 
consistent with those reported for the 102 MHz Geminga detections. We note than even a 
1 pc cm -3 error in our assumed DM would produce only 7 ms of pulse smearing, negligible 
for Geminga's period of 237 ms. 

Each one-hour dedispersed time series was folded using the gamma-ray ephemeris of 
Mattox et al. (1998) (P = 0.23709574610 s, P = 10.974012 x 10~ 15 s s _1 , epoch = 2446600). 
We used a single period to fold each 1 hour data set, but updated the period between scans, 
using the gamma-ray ephemeris and TEMPO ( |l'aylor fc Weisberg 1989] ). We note that 
using a single period for each scan results in only 0.2 ms of pulse smearing. Table 1 lists 
the modified Julian date of the midpoint of each scan, the corresponding period, and the 
total length of the scan. The first 6 scans are 1 hr in length, while the last is 50 minutes, 
yielding a total integration time of 6 hr 50 min. 

Figure 1 shows the folded profiles for scans 1 through 7, in addition to the composite 
profile formed by calculating the phase at the start of each scan, and shifting and summing 
the profiles accordingly. We have verified that our phase- shifting algorithm is correct 
by applying it to the observatory generated 19.2 Hz waveguide switch signal. There is 
structure in some of the profiles of Figure 1. However, no consistent pulse shape and/or 
phase is apparent. Furthermore, coherently summing all profiles dramatically decreases the 
signal-to-noise level of any features. We have also calculated Fourier transforms of all 7 
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dedispersed time series. We have not detected any harmonics of Geminga's 237-ms period. 

Figure 2 shows the folded pulse profiles for our test pulsars, B0329+54 (DM = 26.8 
pc cm -3 ) and B0950+08 (DM = 3.0 pc cm -3 ). These data were processed with identical 
dedispersion and folding routines to the Geminga data. We have subtracted the off-pulse 
mean to calculate F, the phase-averaged pulsed flux density. For PSR B0329+54, F = 1071 
mJy, and for PSR B0950+08, F = 8.9 mJy. We note that we also detected both pulsars 
with high signal-to-noise in Fourier transforms of their time series. 

Because Geminga may be too weak to be detected through its time-averaged flux, we 
searched for Crab-like 'giant' pulses, or individual pulses with amplitudes much greater 
than the mean pulse amplitude ( [Lundgren et al. 1995Q . We search for these aperiodic, 
dispersed pulses by dedispersing each data set over a range of trial DMs and, for each 
dedispersed time series, recording those pulses with amplitudes above some signal-to-noise 
threshold. We enchanced our sensitivity to broadened pulses by repeating this thresholding 
with different levels of time series smoothing. We found no evidence for isolated pulses from 
Geminga. In Figure 3, we show the results of this analysis for Geminga, PSR B0950+08, 
and PSR B0329+54, from which we detect a large number of individual pulses. 

3.1. Bayesian Pulsed Flux Estimation 

In Appendix A, we derive a Bayesian procedure for estimating the pulsed flux density 
given a folded profile with a pulse of unknown amplitude, width, and phase. Assuming the 
additive noise is Gaussian, we are able to calculate a probability density function (PDF) for 
the phase-averaged pulsed flux density independent of pulse width, pulse amplitude, pulse 
phase, off-pulse mean, and off-pulse rms. Because no prior assumptions are necessary to 
calculate an upper limit, this method is an improvement over standard methods used to 
test for pulse shape uniformity and to calculate upper limits. The \ 2 statistic, the method 
most commonly used, returns only the deviation of a binned pulse profile from a uniform 
distribution. More sophisticated tests, such as the Z 2 and H tests (PeJager 1994] ) , developed 
to analyze sparse X-ray and gamma-ray profiles, are less dependent on binning and more 
sensitive to a wide variety of pulse shapes than the \ 2 statistic, but still yield only the 
probability that a distribution differs from uniformity. To calculate an actual upper limit, 
one must still assume a pulse width and phase. Furthermore, these methods do not provide 
a simple mechanism for calculating a PDF for the pulsed flux density. In following papers, 
we will use our method to calculate more accurate X-ray and gamma-ray upper limits for 
known radio pulsars and to search for new high-energy pulsars. The method can also be 
adapted to a wide range of pulse shapes, in addition to the simple square pulse presented 
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in Appendix A. 

We apply the method of Appendix A to calculate PDFs for the pulsed flux density for 
the folded Geminga pulse profiles shown in Figure 1. The resultant PDFs are also shown in 
Figure 1 and are discussed further in Section 3.2.2. 



3.2. Interstellar Scintillations of the Geminga Pulsar 

Since Geminga is a nearby (d ~ 150 pc) pulsar, we expect its flux density to be strongly 
modulated by diffractive interstellar scintillations (DISS). We therefore present a method 
which allows us to calculate a PDF for a source's flux density given several measurements 
of the scintillated flux density. We restrict our analysis to the strong scattering regime (a 
good assumption at low radio frequencies), where rms , the rms phase perturbation due to 
electron-density variations, is much greater than 1. In this case, the amplitude fluctuations 
are 'saturated' because the electric field has been sufficiently randomized to have complex 
Gaussian statistics. 

We describe DISS through a characteristic timescale At d and characteristic bandwidth 
Az/d. These quantities depend upon the scattering measure (SM), the integrated turbulence 
strength (assuming a power spectrum for electron density irregularities) along the line of 
sight to the pulsar. In the strong scattering regime, and assuming a uniform distribution 
of scattering material and a Kolmogorov power spectrum, At d and Ah> d scale with SM, 
observing frequency, distance, and pulsar velocity as flUordes fc Lazio 1991| ; |Gordes & 
Rickett 1998| ) 

At d = lO^V^SM- - 6 ^- 1 s, (1) 

Av d = lO-^V^SM- 1 - 2 ^- 1 kHz (2) 

where v is the radio frequency in GHz, SM is the scattering measure in kpc m _20//3 , V is the 
pulsar's transverse velocity in km s _1 , and D is its distance in kpc. The scaling is such that 
a distant source with high velocity observed at a low frequency has a shorter characteristic 
bandwidth and timescale than a nearby, low velocity source observed at high frequency. 

Once we have characterized the DISS properties of a source (through a measurement 
or an estimation using Eqs. [I] and 0), we calculate the posterior PDF of the intrinsic source 
flux density given several measurements of the scintillated flux density using the method 
described in Appendix B. The source's intrinsic flux density, S, is modulated by interstellar 
scintillations by a factor g, so that the measured flux density, F, is 

F = gS + N, (3) 



- 6- 



where N, the additive noise, may include both radiometer noise and radio frequency 
interference. We note that this additive noise model holds only for the case where the 
time-bandwidth (TB) product of the measurements (ie. integration time T times bandwidth 
B) is large, TB ^> 1. In the strong scattering regime, the scintillation time-bandwidth 
product, AtdAz/d, which is the characteristic size of a "scintle" in the v — t plane, may or 
may not be much larger than TB for the measurement. When the time-bandwidth product 
of the measurement is smaller than that of one DISS 'scintle,' the PDF for g is a one sided 
exponential: 

f 9 (g)=e~ 9 , g>0 (4) 

with CDF 

F g = P{<g} = \- e-°. (5) 



When TB ^> Aid Az/d, the measurements average over many scintles, increasing the 
number of degrees of freedom from 2 (for unquenched DISS) to 2iViss, where 



^iss 



B \ / T 

1 + 0.2 



)( 1 + ' 2 ^)' (6 > 



Then, g is a x 2 random variable with 2iVigg degrees of freedom whose PDF is 

( nNjoo) Njss 

f g (g, Afe) = W™> e -« N °"U(g), (7) 



ISSJ 



where U(g) is the Heaviside function and T is the gamma function. For pulsars at large 
distances or those observed at low frequencies, the quantity iViss will approach infinity and 
f g will tend toward a delta function, S(g — 1). Figure 4 illustrates the dependence of f g on 
AW 



3.2.1. Results for Known Pulsars 

For the known pulsars B0329+54 and B0950+08, Aid and Az^d have been measured 
( |5tinebring et al. 1996| ; Phillips fc Clegg 1992|) . Applying the known frequency scaling of 



these quantities (see Eqs. |T] & 0), we estimate Atd = 142 s, Az/d = 25 kHz for B0329+54 
and At d = 2757 s, Az/ d = 102 MHz for B0950+08 at 317 MHz. Given these values and our 
observing bandwidth {B = 25 MHz) and time (T = 600 s), we may calculate iVigs from 
Eq. |. For B0329+54, iV ISS = 371, TB 3> Az/ d At d , and f g will approach a delta function. 
B0950+08 is much closer to Earth, so its time-bandwidth product Av^At^ > TB and 
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-Wiss = 1-095. Although this pulsar is likely in the transition regime from weak to strong 
scattering, f g will approximate an exponential and will have a similarly long tail. 

Using the method of Appendix B, we calculate the PDFs for the intrinsic source fluxes 
of these pulsars given by our measurements. In Figure 5, we plot fs(S), where S is the 
intrinsic phase-averaged flux density. For B0329+54, the PDF of the intrinsic source flux 
density is strongly peaked at the measured flux density, and follows a roughly Gaussian 
distribution away from the peak. For B0950+08, the PDF for the intrinsic source flux 
density again peaks at the measured flux, but is characterized by a broad exponential 
tail. We can calculate confidence intervals for the flux densities of these pulsars from their 
PDFs. For B0329+54, the 95% confidence interval is 968 mJy < S < 1190 mJy, while for 
B0950+08, the 95% confidence interval is 5.4 mJy < S < 880.1 mJy. We note that these 
pulsars' predicted 317 MHz fluxes from the Princeton Pulsar Catalog fllaylor et al. 199H| ) 
400 MHz fluxes and measured spectral indices ([Lorimer et al. 1995|) are 1890 mJy for PSR 
B0329+54 and 590 mJy for PSR B0950+08. While the predicted flux for PSR B0329+54 
lies outside of our confidence interval, this discrepancy is consistent with the ~ 40% 
modulations due to refractive interstellar scintillation that have been measured for this 
pulsar ( |Stinebring et al. 1996|) . Refractive ISS is probably also important in modulating the 
flux of PSR B0950+08, although its predicted flux does fall within our confidence interval. 



3.2.2. Results for Geming a 

Because the scintillation bandwidth and timescale have not been measured for 
Geminga, we use Eqs. [I] and [2] to predict the characteristic timescale and bandwidth given 
our observing frequency, Geminga's measured distance, and the SM predicted from the TC 
model. For Geminga, we estimate SM fa 10~ 4 ' 4 kpc m -20 / 3 and, from Eq. ^|, Az/d ~ 1.5 
MHz. Given Geminga's transverse velocity of V± ~ 140 km s _1 , Eq. |T] yields At d « 275 
s. Therefore, TB 3> Au d At d , iV ISS « 16, and we expect statistical independence between 
measurements. Applying Eq. [B9] to the distributions for fr(F) (see Figure 1), we calculate 
fs{S) for each of the 7 observations individually. These PDFs are plotted in Figure 6. 



Using Eqs. |B8] and [B£], we may use all 7 Geminga observations to calculate a composite 
PDF, shown in Figure 7. This PDF does peak at a non-zero value of 5*. However, as shown 
in Figure 10 (see Appendix A), the PDFs of simulated noise often peak at non-zero values. 
We note that the rms noise levels in the Geminga pulse profiles are roughly equivalent 
to those of the simulated profiles in this figure. In addition, we expect more structure in 
the Geminga profiles due to the large amount of radio frequency interference at 317 MHz. 
Furthermore, the marginalized PDFs for pulse width and phase for each of the 7 Geminga 
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scans do not peak at a consistent value of these quantities, as expected for a real pulsar. 
We therefore treat our result as an upper limit and calculate a 95% confidence upper limit 
to Geminga's pulsed flux of 3.0 mJy. 

The results so far have assumed that the pulse width is completely unknown. We have 
also calculated upper limits assuming reasonable values for Geminga's pulse width. If we 
assume Geminga's radio pulse has the same width (roughly 180°) as its gamma-ray pulse, 
we calculate a 95% confidence upper limit to the pulsed flux density of 4.0 mJy. However, 
except in the case of the Crab, the radio pulse widths of the known EGRET pulsars are 
much narrower than the gamma-ray pulse widths. For this reason, we also calculate an 
upper limit for a pulse width of 25° (calculated assuming the w oc P -1 / 2 scaling of pulse 
width w with period P given by Biggs (1990)) of 1.6 mJy. 

We note that the distribution of material along the line of sight will most likely not be 
uniform, and Atd and Az/d may be somewhat different from those predicted by Eqs. [1] and 
|2|. We therefore explore what happens in the limiting case where the scintillation timescale 
is longer than the 9-hr total timescale of our observations. For this case, we use Eq. [B^] to 
calculate the PDF for intrinsic source flux density given all 7 observations, again integrating 
over the PDFs of Figure 1. The resultant PDF for this limiting case, also plotted in Figure 
7, is much broader than that derived assuming statistically independent scintillations. The 
95% confidence upper limit on Geminga's pulsed flux for this case is 13.9 mJy. However, 
for the scintillation timescale to be as large as 9 hrs, the scattering measure would have to 
be 10~ 4 times smaller than that estimated using the TC model. Alternatively, a timescale 
this long could also be explained by a thin screen located only 0.021 pc from the observer 
( Cordes fc Rickett 1998| ). Obviously, both of these scenarios are very improbable, and the 
true 95% upper limit is likely very close to 3 mJy. 

4. Discussion and Conclusions 

Comparing our 3 mJy upper limit to Geminga's flux density at 317 MHz with the 
mean flux density of 60 mJy reported by Malofeev & Malov at 102 MHz, we calculate a 
lower limit on the spectral index of Geminga of a ~ 2.7 (F(u) oc v~ a \ comparable to the 
spectral index of the Crab. If the spectral index of Geminga is indeed this high, it would 
be extreme, as the Crab and Geminga pulsars have quite different ages, spin periods, and 
spin-down rates. 

In Figure 8, we plot our upper limit along with previously measured upper limits and 
detections. Unfortunately, while there have been many radio searches for Geminga, there 
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are few published upper limits. This is partly due to Geminga's migrating gamma-ray 
positional error box. For example, Boriakoff et al. (1984) searched for radio pulsars 
in Geminga's error box, determined from COS B data, but did not cover the current, 
well-determined position. Of special interest is the recent non-detection of Geminga at both 
35 and 327 MHz by Ramachandran et al. (1998). While their quoted 327 MHz upper limit 
of 0.3 mJy is much lower than ours, the rms sensitivities of our 7-hr VLA observation and 
their 6-hr Ooty observation are very similiar. For consistency, we have plotted the upper 
limit obtained through applying our method to their data, taking into account scintillation 
and the unknown width and phase of Geminga's pulse. 

The reason for Geminga's radio weakness is still debatable. One explanation is that 
some gamma-ray pulsars simply do not emit at radio energies. Halpern & Ruderman 
(1993) suggest that Geminga's radio emission is quenched by the copious pair production 
in the inner magnetosphere. However, the transient dips in Geminga's soft X-ray profile 
( Halpern fc Wang 19971 ) suggest that the processes responsible for supplying e ± to the 
inner magnetosphere may be variable. Because the radio silence will be occasionally broken 
as this e ± quenching plasma clears away, we may expect Geminga to be a transient radio 
emitter. If Geminga continues to be the only high-energy pulsar found to exhibit this 
transient phenomenon, studying its X-ray properties may be useful for determining if and 
why radio emission is suppressed. 

A more likely explanation for Geminga's radio weakness is simply that there is some 
misalignment between the gamma-ray and radio pulsar beams. Geminga's pulse shape in 
gamma-rays (like that of the Crab and Vela pulsars) is broad and double-peaked, suggesting 
an origin in the "outer gaps" of the pulsar magnetosphere near the light cylinder. Because 
radio emission is expected to be associated with the open field line region centered on the 
magnetic poles, it is possible that our line-of-sight intersects the gamma-ray beam only. 
Romani & Yadigaroglu (1995), using the outer gap model for pulsar gamma-ray emission, 
have modeled the orientation and size of the radio and gamma-ray pulsar beams given radio 
and gamma-ray data on known pulsars. They find that, because the gamma-ray beam 
is much wider than the radio beam, 45% of young pulsars will be detected only at high 
energies, while only 19% of young pulsars will be detected at both radio and gamma-ray 
energies. Scaling from the 5 EGRET detections of radio pulsars, Romani & Yadigaroglu 
expect 12 nonradio pulsars to be visible in gamma-rays at flux levels comparable to the 
radio selected objects. This implies that most of the ~ 20 unidentified Galactic EGRET 
sources are young radio-quiet or radio-weak pulsars like Geminga. This is consistent with 
other studies which find that the properties, such as flux variability ( |McLaughlin et al" 



1996| ) and spectral index ( [Merck et al. 199"6| ), of many unidentified sources are similar to 



those of the known EGRET pulsars. It has also been shown ([Yadigaroglu fc Romani 19971 ) 



-10- 



that many of the unidentified sources are associated with regions of massive star formation 
and death, which are expected to be breeding grounds for pulsars. 
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A. Bayesian Pulsed Flux Estimation 

Pulsar flux densities are usually reported as phase-averaged quantities. Therefore, to 
calculate a pulsed flux upper limit, one must assume a pulse width, pulse phase, off-pulse 
mean, and off-pulse rms. To avoid trial and error estimates, we have developed a method 
to calculate a PDF for the pulsed flux density independent of these quantities. 

We assume an M-bin folded pulse profile, where the amplitude of each bin is described 
by di (i e 1, M). Our model for the expected amplitude within each bin, is parametrized 
by 

w - width of square pulse 

m - first bin of pulse 

F - phase-averaged pulsed flux 

/i - off-pulse mean 

cr - off-pulse rms 

Given this model, we may express di as 

- f FMI w + Ni i = m, m+w— 1 . . , . 

di = < AT . (Al) 

[ Ni otherwise 

where Ni is the amplitude of the noise in bin i. The likelihood function of the data is given 
by the joint probability density function (PDF) of the noise samples. Assuming Gaussian 
noise with mean p, and rms a, and statistical independence of the bins, we may write this as 



£ = (2na 2 )- M/2 



cxp 



1 M 



2^ 2 l=1 



(A2) 



To find the best model parameters, we use Bayes' Theorem to write the posterior 
probability p(F, w, m, a\d, I) of a model parametrized by F, w, m,p:,Sza given the data 
d and any prior information I as 

p(d\F,w,m,n,a,I)p(F,w,m,/i,a\T) 

p{F, w, m, y., a d, I) = — — . (A3) 

p[d\T) 



The factor p(d\F, m, w, fi, a, I) is simply the likelihood function given in Eq. |A2| . Because we 
expect all parameter values to be equally likely a priori, we assume flat priors for F, w, m, 
& ji (i.e. p(F\T) = constant) and a prior oc 1/a (i.e. p(log(cr)|I) = constant). We 

choose a prior that is uniform in logo" to express our ignorance of the scale of the noise 
variation. Then, we can express the posterior probability as simply 

p(F, m, w, fi, a\d, I) oc — . (A4) 

a 



-12- 



Since we would like to calculate the PDF of the phase- averaged flux F independent 
of the other parameters, we must marginalize the posterior probability over the "nuisance 
parameters" w,m,fi, & a. We write the PDF for F as 

fp{F) oc J — J dm J dfx J da p(F, w, m, /i, a\d, I). (A5) 



Substituting Eq. [A^ into Eq. |A5| , we find that the integrals over fx and a may be done 
analytically. After some algebra, we find 

, fdwr, / F 2 M 2 2FMD P (FM-A) 2 \ Mf)/2 

f F F oc / — dm [D 2 + p -- y - rp^M , A6 

J w J \ w w M J 

where constant factors have been eliminated and 

M M m+w—1 

D x = Y i d i , D 2 = "£dl & D p = d i- (A7) 

i=l i=l i=m 

In the case that the width and/or phase of the pulse is known, the marginalization over 
m and/or w in Eq. M can be removed to further constrain the PDF for F. 



To test our method, we have created many simulated profiles consisting of a pulse 
superimposed on Gaussian noise, and have confirmed that the PDF of Eq. |A6| does 
maximize at the correct value of F. Some example simulation results for square wave pulses 
are shown in Figure 9. Figure 10 shows the results of applying our algorithm to profiles 
consisting only of Gaussian random noise. 



B. Bayesian Flux Estimation of a Scintillating Source 

We assume that a source's intrinsic flux density is modulated by interstellar 
scintillations by a factor g, so that the measured flux density is 

F = gS + N, (Bl) 

where is the additive noise and S is the intrinsic flux density. To calculate a PDF for S 
we assume we have m measurements of F and m separate estimates of N, corresponding to 
making on and off-pulse flux measurements. We assume that the noise PDF is Gaussian (a 
good assumption for radio observations with large time-bandwidth product), and that we 
have estimated ii, the average value of the off-pulse mean, and a, the standard deviation 
of this quantity, for each of m profiles. The PDF of the phase-averaged noise, f^(N), 
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is a Gaussian with mean \i and standard deviation a, and the PDF for an individual 
measurement of F is 

f F (F) = J dNf N (N)f 9 (^-^) , (B2) 

where f g is the PDF of g. The likelihood function is a combination of all the data using 
PDF's of this form. So, for m statistically independent observations, the likelihood function 
for S is 

m 

£(S) = I] fr{Fj) (B3) 

3=1 

We cannot use this scheme for all cases because generally we need to use the joint 
PDF for the scintillation gain g for the different measurements. Because there is no closed 
expression for this, we consider only two limiting cases. 



B.l. Case I: Statistically Independent ISS 

The first case we consider is when the interstellar scintillations are statistically 
independent between measurements. For this case, we may write 

WS) = 11/ dNf N (N)f g (-^— ) ■ (B4) 

where /jv(AT) is a Gaussian described by fij and Oj. 

For Atd "C T (i.e. the ISS time scale Atd is much less than the averaging time T) 
and/or Az/ d -C B (the ISS bandwidth is much less than the receiver bandwidth), the ISS is 
quenched and has a PDF that is related to a x 2 PDF with the number of degrees of freedom 
given by 2iViss, where Niss is the number of scintles averaged over. In this case, statistical 
independence between observations is a good assumption. 



B.2. Case II: Perfectly Correlated Scintillations 

The second case we consider is when the time scale for interstellar scintillations is much 
longer than the total data acquisition time (e.g. mT for contiguous measurements). In 
this case, the scintillation gain g is identical for all measurements, and has some unknown 
value with PDF described by Eq. In this case, individual measurements are statistically 
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independent for the noise but completely identical with respect to g. The likelihood function 
is therefore simply the product of the noise PDFs, and can be expressed as 

m m 

£s(S) = II fpiFj) = II In{F 3 - S). (B5) 
3=1 3=1 

We can use this equation to get the likelihood function for the product 5" = gS and then 
integrate over f g (g) to get the likelihood for S: 

C S (S) = Jdgf g (g)C s ,(gS) (B6) 

m 

where U(S') = flMFj-S'). (B7) 



Once we have chosen the scintillation regime of interest, and calculated Cs(S) with 
Eq. [B4] or Eq. |B7| , we may calculate the PDF of S as the normalized likelihood function 

fs{S) = JS°dSC 8 (S)- (B8) 

When no source contributes to the 'on-source' measurement, we expect fs(S) to 
maximize at S = and have a width determined by a and the number of measurements. 
An upper bound on S can be calculated by choosing a probability level 1 — e and calculating 
the value of S such that the area of fs{S) above that value is e. If fs{S) maximizes at a 
non-zero value, then a confidence interval for S can be similarly calculated. 

When, as in the calculation of an upper limit using the method of Appendix A, we 
do not have a single measured value for Fj, but instead a PDF, we must replace Fj by an 
integral over its PDF. In this case, Eq. |B3] becomes 

m r r /F — N\ 

£s(S) = n J fp{Fj)dFj j dNf N (N)f g [- L ^—) ■ (B9) 
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Table 1. 



Scan number 


MJD 


Period (s) 


Time (s) 


1 


50847.08182 


0.237113398987 


3600 


2 


50847.13715 


0.237113501528 


3600 


3 


50847.18646 


0.237113605070 


3600 


4 


50847.23252 


0.237113704454 


3600 


5 


50847.27766 


0.237113797450 


3600 


6 


50847.34061 


0.237113908108 


3600 


7 


50847.38646 


0.237113967595 


3000 



0.25 0.5 0.75 1 

Pulse Phase (cycles) 



2 4 6 8 10 

F (mJy) 



Fig. 1. — Left: The upper 7 plots show the 64 bin folded pulse profiles for the individual 
data sets. The y-axis has units of mJys (one tick mark = 10 mJy) and have been normalized 
to zero mean. The profiles have been shifted so that the first bin of each profile corresponds 
to the same pulse phase. The bottom plot shows the profiles from all 7 individual phase- 
referenced scans summed. Right: The upper 7 plots show /f(F) for the individual Geminga 
scans. The lower plot shows /f(-F) for the summed profile. 
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Fig. 2.— Folded pulse profiles for PSR B0329+54 (left) and PSR B0950+08 (right) are 
shown. The off-pulse mean has been subtracted. The phase-averaged pulsed flux density is 
1071 mJy for B0329+54 and 8.9 mJy for B0950+08. 



-20- 




Fig. 3. — Upper: For PSR B0329+54, all individual pulses with amplitudes above 9er have 
been plotted vs. DM and arrival time. This pulsar, at a DM of 26.8 pc cm -3 , is clearly 
detectable through its individual pulses. Lower: The number of isolated pulses above a 4 
a threshold is plotted vs. DM for (from left to right) B0329+54, B0950+08, and Geminga. 
While a peak at the DM of PSR B0329+54 is again obvious, there is no evidence of individual 
pulses from PSR B0950+08 or Geminga, both with DM m 3 pc cm~ 3 . The individual pulses 
detected in these directions are most likely due to the considerable interference at 317 MHz. 
We note that the baseline of the histograms is not zero, and the distributions of pulses for 
PSR B0950+08 and Geminga are consistent with a uniform distribution in DM. We also 
note that, while the number of pulses above threshold is 25 times greater for Geminga than 
for B0950+08, the total time of the Geminga observations is greater by a factor of 40, and 
so the rate of pulse detection is lower for Geminga. 
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Fig. 4. — The PDF of g is plotted for different values of iV ISS , defined in Eq. 6. The lines 
shown are for (from thinnest to thickest) iV^s = 1,1.5,3,6,10, and are tagged with their 
iViss values. 
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S (mJy) 



Fig. 5. — The PDFs for phase-averaged intrinsic source flux density S, given the measured 
phase-averaged fluxes (from Figure 2), for our two test pulsars are shown. For PSR 
B0329+54, Niss = 371, and for PSR B0950+08, N lss = 1.095. 
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Fig. 6. — The upper 7 plots show fs(S) for the individual Geminga scans. The lower plot 
shows fs(S) for the summed profile. 



1 2 3 4 5 



5 10 15 20 



S (mJy) 

Fig. 7. — The left plot shows the composite PDF for source flux density S calculated from all 
7 Geminga scans, assuming statistical independence among measurements. The right plot 
shows the composite PDF for S calculated in the limiting case where the timescale for ISS 
is longer than the total time of our observation. The dotted lines mark the 95% confidence 
upper limits. 
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Fig. 8. — This plot shows all published Geminga upper limits and detections, along with 
error bars (when available). The three 102.5 MHz detections are offset slightly from each 
other for clarity. The slope of the solid line (2.65) represents our calculated lower limit on 
Geminga's spectral index. 
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Fig. 9. — The top plots show two simulated pulse profiles with M = 64, // = 0, and a = 10. 
For profile (1), w — 32 and F = 30, and for profile (2), w — 6 and F = 5.625. Applying the 
Bayesian pulse detection scheme to these profiles results in the PDFs shown in the bottom 
plots. Note that these PDFs maximize at F, with scatter determined by the Gaussian noise 
in the original profiles. 
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Fig. 10. — This plot show the PDFs resulting from application of the Bayesian pulse detection 
scheme to 64 bin profiles of only Gaussian noise with /i — and a = 10. The results of 100 
independent trials are plotted. Some PDFs maximize at non-zero values, all less than the 
rms noise level. 



